Chapter 4 Community composition

load("data/data.Rdata")

4.1 Taxonomy overview

4.1.1 Stacked barplot

# Merge data frames based on sample
merged_data<-genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
  filter(count > 0) #filter 0 counts

# Create an interaction variable for time_point and sample
merged_data$interaction_var <- interaction(merged_data$sample, merged_data$time_point)

ggplot(merged_data, aes(x=interaction_var,y=count, fill=phylum, group=phylum)) + #grouping enables keeping the same sorting of taxonomic units
    geom_bar(stat="identity", colour="white", linewidth=0.1) + #plot stacked bars with white borders
    scale_fill_manual(values=phylum_colors) +
    facet_nested(. ~ type,  scales="free") + #facet per day and treatment
    guides(fill = guide_legend(ncol = 1)) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
          axis.title.x = element_blank(),
          panel.background = element_blank(),
          panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black")) +
   labs(fill="Phylum",y = "Relative abundance",x="Samples")+
  scale_x_discrete(labels = function(x) gsub(".*\\.", "", x)) +
  facet_wrap(~type, scales = "free", labeller = as_labeller(function(label) gsub(".*\\.", "", label))) #only show time_point label

4.1.2 Phylum relative abundances

phylum_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  left_join(genome_metadata, by = join_by(genome == genome)) %>%
  group_by(sample,phylum) %>%
  summarise(relabun=sum(count))

phylum_summary %>%
    group_by(phylum) %>%
    summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
    arrange(-mean) %>%
    tt()
tinytable_rnwtr57oxrkslwyxnrzc
phylum mean sd
p__Bacteroidota 0.380806816 0.202118047
p__Bacillota_A 0.253813952 0.163904042
p__Bacillota 0.118079729 0.146679823
p__Pseudomonadota 0.096901291 0.160581847
p__Campylobacterota 0.053050348 0.092998130
p__Verrucomicrobiota 0.029867960 0.073128363
p__Desulfobacterota 0.023580475 0.036482484
p__Chlamydiota 0.010572704 0.059690269
p__Fusobacteriota 0.010482132 0.028311738
p__Cyanobacteriota 0.009206465 0.016492133
p__Bacillota_C 0.004713164 0.006645703
p__Spirochaetota 0.004009680 0.012308028
p__Bacillota_B 0.002465465 0.004927667
p__Actinomycetota 0.001235741 0.006346852
p__Elusimicrobiota 0.001214079 0.006501752
phylum_arrange <- phylum_summary %>%
    group_by(phylum) %>%
    summarise(mean=mean(relabun)) %>%
    arrange(-mean) %>%
    select(phylum) %>%
    pull()

phylum_summary %>%
    filter(phylum %in% phylum_arrange) %>%
    mutate(phylum=factor(phylum,levels=rev(phylum_arrange))) %>%
    ggplot(aes(x=relabun, y=phylum, group=phylum, color=phylum)) +
        scale_color_manual(values=phylum_colors[rev(phylum_arrange)]) +
        geom_jitter(alpha=0.5) + 
        theme_minimal() + 
        theme(legend.position="none") +
        labs(y="Phylum",x="Relative abundance")

4.2 Taxonomy boxplot

4.2.1 Family

family_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  group_by(sample,family) %>%
  summarise(relabun=sum(count))

family_summary %>%
    group_by(family) %>%
    summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
    arrange(-mean) %>%
    tt()
tinytable_2btwr2xncpn8tu5tuia9
family mean sd
f__Bacteroidaceae 2.215146e-01 0.1392048812
f__Lachnospiraceae 1.406192e-01 0.1085410432
f__Tannerellaceae 1.028077e-01 0.0799387840
f__Helicobacteraceae 5.258990e-02 0.0925721039
f__Mycoplasmoidaceae 3.694491e-02 0.0756423660
f__Erysipelotrichaceae 3.502391e-02 0.0451216980
f__UBA3700 3.418595e-02 0.0557485850
f__Marinifilaceae 2.774048e-02 0.0271999075
f__ 2.772765e-02 0.0838376915
f__DTU072 2.700317e-02 0.0529029477
f__Rikenellaceae 2.696900e-02 0.0464477971
f__Enterobacteriaceae 2.691677e-02 0.0918386902
f__Coprobacillaceae 2.551773e-02 0.0892179821
f__Desulfovibrionaceae 2.358047e-02 0.0364824837
f__Ruminococcaceae 1.859938e-02 0.0429483670
f__LL51 1.759437e-02 0.0687507595
f__Rhizobiaceae 1.596807e-02 0.0770633474
f__UBA3830 1.578817e-02 0.0454016753
f__Akkermansiaceae 1.227359e-02 0.0316838689
f__Chlamydiaceae 1.057270e-02 0.0596902690
f__Fusobacteriaceae 1.048213e-02 0.0283117384
f__CAG-239 9.141331e-03 0.0150992357
f__Enterococcaceae 8.420523e-03 0.0466580889
f__Gastranaerophilaceae 7.728770e-03 0.0144404459
f__Oscillospiraceae 6.643218e-03 0.0078105983
f__UBA1997 6.378408e-03 0.0309668631
f__Streptococcaceae 6.366441e-03 0.0342174908
f__UBA1242 4.673359e-03 0.0153730061
f__Brevinemataceae 4.009680e-03 0.0123080282
f__Acutalibacteraceae 3.374550e-03 0.0109560112
f__UBA660 3.196511e-03 0.0139558713
f__Clostridiaceae 2.842205e-03 0.0171339212
f__RUG11792 2.817729e-03 0.0251127757
f__Peptococcaceae 2.465465e-03 0.0049276669
f__MGBC116941 2.160169e-03 0.0093827415
f__Acidaminococcaceae 1.943804e-03 0.0050321931
f__CAG-508 1.807978e-03 0.0064175379
f__Anaerovoracaceae 1.569531e-03 0.0036471440
f__Moraxellaceae 1.485415e-03 0.0097446552
f__RUG14156 1.477695e-03 0.0045847831
f__Staphylococcaceae 1.361709e-03 0.0050863070
f__Elusimicrobiaceae 1.214079e-03 0.0065017516
f__CAG-288 9.490865e-04 0.0060198775
f__Anaerotignaceae 8.989745e-04 0.0040469504
f__CALVMC01 7.516846e-04 0.0043609744
f__Eggerthellaceae 6.407882e-04 0.0021266467
f__Massilibacillaceae 6.235073e-04 0.0016350480
f__Mycobacteriaceae 5.949530e-04 0.0060400054
f__UBA1820 4.736030e-04 0.0013057353
f__Arcobacteraceae 4.604457e-04 0.0050324221
f__CAG-274 4.519746e-04 0.0022028477
f__Burkholderiaceae_C 3.699430e-04 0.0048092594
f__Muribaculaceae 3.617894e-04 0.0009776409
f__UBA932 3.419549e-04 0.0011583178
f__Hepatoplasmataceae 2.989107e-04 0.0038858387
f__Rhodobacteraceae 2.959092e-04 0.0038468195
f__Weeksellaceae 2.771627e-04 0.0031476408
f__Eubacteriaceae 1.646822e-04 0.0006729067
f__Sphingobacteriaceae 1.505775e-04 0.0012460017
f__Devosiaceae 1.489995e-04 0.0015094318
f__Pumilibacteraceae 1.277418e-04 0.0007646754
f__WRAU01 9.603359e-05 0.0012484367
f__Peptostreptococcaceae 2.287338e-05 0.0002973539
family_arrange <- family_summary %>%
    group_by(family) %>%
    summarise(mean=sum(relabun)) %>%
    arrange(-mean) %>%
    select(family) %>%
    pull()

family_summary %>%
    left_join(genome_metadata %>% select(family,phylum) %>% unique(),by=join_by(family==family)) %>%
    left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
    filter(family %in% family_arrange[1:20]) %>%
    mutate(family=factor(family,levels=rev(family_arrange[1:20]))) %>%
    filter(relabun > 0) %>%
    ggplot(aes(x=relabun, y=family, group=family, color=phylum)) +
        scale_color_manual(values=phylum_colors[-8]) +
        geom_jitter(alpha=0.5) + 
        facet_grid(.~type)+
        theme_minimal() + 
        labs(y="Family", x="Relative abundance", color="Phylum")

4.2.2 Genus

genus_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
  left_join(genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  group_by(sample,genus) %>%
  summarise(relabun=sum(count)) %>%
  filter(genus != "g__")

genus_summary %>%
    group_by(genus) %>%
    summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
    arrange(-mean) %>%
    tt()
tinytable_hy8zrnj1dtq6zlpooa2i
genus mean sd
g__Bacteroides 1.352836e-01 0.0923254347
g__Parabacteroides 9.681305e-02 0.0801790006
g__Phocaeicola 6.935126e-02 0.0794715561
g__Mycoplasmoides 3.051289e-02 0.0753256189
g__Helicobacter_J 3.009004e-02 0.0594808174
g__Odoribacter 2.552307e-02 0.0267570773
g__Roseburia 2.384718e-02 0.0559439989
g__NHYM01 2.249986e-02 0.0796951054
g__Alistipes 2.208351e-02 0.0283591693
g__Coprobacillus 2.013859e-02 0.0878851683
g__Agrobacterium 1.596807e-02 0.0770633474
g__Akkermansia 1.227359e-02 0.0316838689
g__Fusobacterium_A 1.038903e-02 0.0283170720
g__Kineothrix 8.801109e-03 0.0409037874
g__Proteus 8.657984e-03 0.0681831315
g__Dielma 8.578108e-03 0.0090836958
g__CAG-95 8.109350e-03 0.0205019603
g__JAAYNV01 7.994602e-03 0.0195404994
g__UBA866 7.260325e-03 0.0293499698
g__Desulfovibrio 7.018155e-03 0.0211479806
g__Enterococcus 7.001300e-03 0.0456600395
g__Ureaplasma 6.432013e-03 0.0138012183
g__Lactococcus 6.366441e-03 0.0342174908
g__Lacrimispora 6.064761e-03 0.0098190433
g__Parabacteroides_B 5.994680e-03 0.0100174097
g__CALXRO01 5.765728e-03 0.0308524398
g__Citrobacter 5.717035e-03 0.0334547978
g__NSJ-61 5.560111e-03 0.0198839325
g__Breznakia 5.504528e-03 0.0237016548
g__Clostridium_AQ 5.326190e-03 0.0121695008
g__Salmonella 5.133583e-03 0.0184556360
g__Bilophila 4.983840e-03 0.0088456017
g__Hungatella_A 4.811802e-03 0.0095549574
g__MGBC136627 4.364796e-03 0.0163003355
g__Escherichia 4.188365e-03 0.0266100578
g__UMGS1251 4.159842e-03 0.0072716746
g__Hungatella 4.116388e-03 0.0190788279
g__Clostridium_Q 4.011996e-03 0.0052127439
g__Brevinema 4.009680e-03 0.0123080282
g__Thomasclavelia 3.902580e-03 0.0109042017
g__Scatousia 3.649941e-03 0.0102727788
g__Enterocloster 3.620384e-03 0.0047339282
g__Mailhella 3.612079e-03 0.0102470769
g__Copromonas 3.599368e-03 0.0050180473
g__Ventrimonas 3.513355e-03 0.0071233140
g__Caccovivens 3.343316e-03 0.0122194941
g__Lawsonia 3.289015e-03 0.0117180962
g__Fournierella 3.217226e-03 0.0062309921
g__Limenecus 3.163279e-03 0.0065731820
g__MGBC133411 3.042418e-03 0.0074576633
g__Mucinivorans 2.900095e-03 0.0373193950
g__Sarcina 2.842205e-03 0.0171339212
g__Acetatifactor 2.738138e-03 0.0058449483
g__Eisenbergiella 2.697412e-03 0.0068331645
g__Bacteroides_G 2.682722e-03 0.0346150378
g__CAJLXD01 2.633818e-03 0.0087495783
g__Blautia 2.558708e-03 0.0061407647
g__C-19 2.275515e-03 0.0048691968
g__Butyricimonas 2.217402e-03 0.0050081070
g__Velocimicrobium 2.201224e-03 0.0066764022
g__CAZU01 2.196385e-03 0.0065944622
g__MGBC116941 2.160169e-03 0.0093827415
g__Intestinimonas 2.081898e-03 0.0039382613
g__Negativibacillus 2.069077e-03 0.0055137501
g__Rikenella 1.985389e-03 0.0037067264
g__Phascolarctobacterium 1.943804e-03 0.0050321931
g__RGIG6463 1.789843e-03 0.0039682804
g__JALFVM01 1.742360e-03 0.0038623852
g__Oscillibacter 1.510459e-03 0.0024992661
g__Pseudoflavonifractor 1.502976e-03 0.0027706264
g__Acinetobacter 1.485415e-03 0.0097446552
g__Citrobacter_A 1.392923e-03 0.0060348874
g__Staphylococcus 1.361709e-03 0.0050863070
g__RGIG4733 1.303790e-03 0.0040480788
g__UBA1436 1.214079e-03 0.0065017516
g__Lachnotalea 1.208197e-03 0.0049169585
g__Ruthenibacterium 1.202021e-03 0.0053632179
g__14-2 1.184643e-03 0.0096299356
g__Beduini 1.173950e-03 0.0025142208
g__Scatocola 1.121193e-03 0.0044975564
g__Faecisoma 1.085585e-03 0.0055596648
g__Enterococcus_A 1.083991e-03 0.0099150522
g__Faecimonas 9.880526e-04 0.0083079244
g__RGIG9287 9.638826e-04 0.0093228400
g__CAG-345 9.490865e-04 0.0060198775
g__Blautia_A 9.208028e-04 0.0029082304
g__RGIG1896 8.343610e-04 0.0062772562
g__CAG-269 7.982448e-04 0.0047176841
g__Marseille-P3106 7.938268e-04 0.0017618833
g__WRHT01 6.429563e-04 0.0026979819
g__Eggerthella 6.407882e-04 0.0021266467
g__IOR16 6.398942e-04 0.0020651222
g__Anaerotruncus 6.265445e-04 0.0016948517
g__RUG14156 6.151780e-04 0.0022147136
g__CHH4-2 6.145042e-04 0.0019997615
g__Corynebacterium 5.949530e-04 0.0060400054
g__Serratia_A 5.860616e-04 0.0076188009
g__CAG-56 4.915432e-04 0.0016341973
g__Merdimorpha 4.736030e-04 0.0013057353
g__MGBC140009 4.679334e-04 0.0024067320
g__CALURL01 4.635287e-04 0.0016737486
g__Aliarcobacter 4.604457e-04 0.0050324221
g__Scatenecus 4.520215e-04 0.0019760876
g__RGIG8482 4.399064e-04 0.0029753981
g__Enterobacter 4.073437e-04 0.0041317730
g__Klebsiella 4.054439e-04 0.0048910857
g__Caccenecus 3.941198e-04 0.0017802371
g__Alcaligenes 3.699430e-04 0.0048092594
g__Plesiomonas 3.633249e-04 0.0027105056
g__HGM05232 3.617894e-04 0.0009776409
g__JAHHSE01 3.592120e-04 0.0014900895
g__Egerieousia 3.419549e-04 0.0011583178
g__Emergencia 3.413630e-04 0.0017450341
g__Enterococcus_B 3.352316e-04 0.0022266910
g__Stoquefichus 3.026072e-04 0.0020503969
g__Hepatoplasma 2.989107e-04 0.0038858387
g__Paracoccus 2.959092e-04 0.0038468195
g__Moheibacter 2.771627e-04 0.0031476408
g__Scatomorpha 2.641015e-04 0.0010184339
g__UBA7185 2.434328e-04 0.0014558191
g__Eubacterium 1.646822e-04 0.0006729067
g__Sphingobacterium 1.505775e-04 0.0012460017
g__Devosia 1.489995e-04 0.0015094318
g__Anaerosporobacter 1.454112e-04 0.0012747262
g__Caccomorpha 1.383122e-04 0.0010540604
g__UBA2658 1.307451e-04 0.0007204965
g__Protoclostridium 1.277418e-04 0.0007646754
g__Angelakisella 1.268479e-04 0.0009221501
g__Cetobacterium_A 9.310217e-05 0.0008718575
g__Rahnella 6.470705e-05 0.0008411917
g__Peptostreptococcus 2.287338e-05 0.0002973539
genus_arrange <- genus_summary %>%
    group_by(genus) %>%
    summarise(mean=sum(relabun)) %>%
    filter(genus != "g__")%>%
    arrange(-mean) %>%
    select(genus) %>%
    mutate(genus= sub("^g__", "", genus)) %>%
    pull()

genus_summary %>%
    left_join(genome_metadata %>% select(genus,phylum) %>% unique(),by=join_by(genus==genus)) %>%
    left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
    mutate(genus= sub("^g__", "", genus)) %>%
    filter(genus %in% genus_arrange[1:20]) %>%
    mutate(genus=factor(genus,levels=rev(genus_arrange[1:20]))) %>%
    filter(relabun > 0) %>%
    ggplot(aes(x=relabun, y=genus, group=genus, color=phylum)) +
        scale_color_manual(values=phylum_colors[-c(3,4,6,8)]) +
        geom_jitter(alpha=0.5) + 
        facet_grid(.~type)+
        theme_minimal() + 
        labs(y="Family", x="Relative abundance", color="Phylum")

4.3 Alpha diversity

# Calculate Hill numbers
richness <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 0) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(richness = 1) %>%
  rownames_to_column(var = "sample")

neutral <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(neutral = 1) %>%
  rownames_to_column(var = "sample")

phylogenetic <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, tree = genome_tree) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(phylogenetic = 1) %>%
  rownames_to_column(var = "sample")

# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

functional <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, dist = dist) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(functional = 1) %>%
  rownames_to_column(var = "sample") %>%
  mutate(functional = if_else(is.nan(functional), 1, functional))

# Merge all metrics
alpha_div <- richness %>%
  full_join(neutral, by = join_by(sample == sample)) %>%
  full_join(phylogenetic, by = join_by(sample == sample)) %>%
  full_join(functional, by = join_by(sample == sample))

4.3.1 Wild samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="0_Wild") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b", "#6b7398")) +
      scale_fill_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b50", "#6b739850")) +
      facet_wrap(. ~ metric ) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.2 Acclimation samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="1_Acclimation") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b", "#6b7398")) +
      scale_fill_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b50", "#6b739850")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.3 Antibiotics samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="2_Antibiotics") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b", "#6b7398")) +
      scale_fill_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b50", "#6b739850")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.4 Transplant_1 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="3_Transplant1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b183","#d57d2c", "#4477AA")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b18350","#d57d2c50", "#4477AA50")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.5 Transplant_2 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="4_Transplant2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b183","#d57d2c", "#4477AA")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b18350","#d57d2c50", "#4477AA50")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.6 Post-Transplant_1 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="5_Post-FMT1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b183","#d57d2c", "#4477AA")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b18350","#d57d2c50", "#4477AA50")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.7 Post-Transplant_2 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b183","#d57d2c", "#4477AA")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b18350","#d57d2c50", "#4477AA50")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.4 Beta diversity

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, tree = genome_tree)

beta_q1f <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, dist = dist)

4.5 Permanovas

4.5.0.1 Load required data

meta <- column_to_rownames(sample_metadata, "Tube_code")

4.5.1 1. Are the wild populations similar?

4.5.1.1 Wild: P.muralis vs P.liolepis

wild <- meta %>%
  filter(time_point == "0_Wild")

# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

wild.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(wild))]
identical(sort(colnames(wild.counts)), sort(as.character(rownames(wild))))

wild_nmds <- sample_metadata %>%
  filter(time_point == "0_Wild")

4.5.1.2 Number of samples used

[1] 28
beta_div_richness_wild<-hillpair(data=wild.counts, q=0)
beta_div_neutral_wild<-hillpair(data=wild.counts, q=1)
beta_div_phylo_wild<-hillpair(data=wild.counts, q=1, tree=genome_tree)
beta_div_func_wild<-hillpair(data=wild.counts, q=1, dist=dist)

4.5.1.3 Richness

Df SumOfSqs R2 F Pr(>F)
species 1 1.631953 0.207198 6.795072 0.001
Residual 26 6.244347 0.792802 NA NA
Total 27 7.876300 1.000000 NA NA

4.5.1.4 Neutral

Df SumOfSqs R2 F Pr(>F)
species 1 2.018797 0.2566342 8.976049 0.001
Residual 26 5.847641 0.7433658 NA NA
Total 27 7.866438 1.0000000 NA NA

4.5.1.5 Phylogenetic

Df SumOfSqs R2 F Pr(>F)
species 1 0.3786052 0.2108638 6.947419 0.001
Residual 26 1.4168908 0.7891362 NA NA
Total 27 1.7954960 1.0000000 NA NA

4.5.1.6 Functional

Df SumOfSqs R2 F Pr(>F)
species 1 0.0787916 0.1594096 4.930642 0.064
Residual 26 0.4154800 0.8405904 NA NA
Total 27 0.4942716 1.0000000 NA NA
Run 0 stress 0.1097679 
Run 1 stress 0.09968339 
... New best solution
... Procrustes: rmse 0.03851713  max resid 0.1653681 
Run 2 stress 0.1270767 
Run 3 stress 0.1631598 
Run 4 stress 0.09968327 
... New best solution
... Procrustes: rmse 0.0001782053  max resid 0.0006750373 
... Similar to previous best
Run 5 stress 0.09968332 
... Procrustes: rmse 0.0001274813  max resid 0.0004797134 
... Similar to previous best
Run 6 stress 0.1128823 
Run 7 stress 0.1262094 
Run 8 stress 0.09968334 
... Procrustes: rmse 0.0001049072  max resid 0.0003900285 
... Similar to previous best
Run 9 stress 0.1150472 
Run 10 stress 0.09968337 
... Procrustes: rmse 0.0001678799  max resid 0.0006352602 
... Similar to previous best
Run 11 stress 0.1097931 
Run 12 stress 0.09968329 
... Procrustes: rmse 8.479073e-05  max resid 0.0003191259 
... Similar to previous best
Run 13 stress 0.09968327 
... Procrustes: rmse 7.892079e-06  max resid 2.618268e-05 
... Similar to previous best
Run 14 stress 0.1262094 
Run 15 stress 0.1128825 
Run 16 stress 0.09968329 
... Procrustes: rmse 5.286762e-05  max resid 0.000195175 
... Similar to previous best
Run 17 stress 0.1262098 
Run 18 stress 0.1097933 
Run 19 stress 0.1023646 
Run 20 stress 0.1128821 
*** Best solution repeated 7 times
Run 0 stress 0.009648577 
Run 1 stress 0.01181313 
Run 2 stress 0.01419841 
Run 3 stress 0.01079754 
Run 4 stress 0.02069535 
Run 5 stress 0.01235487 
Run 6 stress 0.0123458 
Run 7 stress 0.009487685 
... New best solution
... Procrustes: rmse 0.01679083  max resid 0.0627009 
Run 8 stress 0.01446298 
Run 9 stress 0.01205477 
Run 10 stress 0.01187749 
Run 11 stress 0.01663114 
Run 12 stress 0.01434102 
Run 13 stress 0.01419545 
Run 14 stress 0.01377079 
Run 15 stress 0.0162717 
Run 16 stress 0.01439811 
Run 17 stress 0.01616874 
Run 18 stress 0.01636192 
Run 19 stress 0.01205573 
Run 20 stress 0.01431334 
Run 21 stress 0.02163577 
Run 22 stress 0.01193592 
Run 23 stress 0.01777682 
Run 24 stress 0.01178341 
Run 25 stress 0.01495281 
Run 26 stress 0.01817262 
Run 27 stress 0.01022689 
Run 28 stress 0.01629507 
Run 29 stress 0.009872308 
... Procrustes: rmse 0.02337549  max resid 0.09257717 
Run 30 stress 0.009968707 
... Procrustes: rmse 0.02613689  max resid 0.100782 
Run 31 stress 0.01174334 
Run 32 stress 0.009521427 
... Procrustes: rmse 0.01195821  max resid 0.04404765 
Run 33 stress 0.01411711 
Run 34 stress 0.01252358 
Run 35 stress 0.01000363 
Run 36 stress 0.01131051 
Run 37 stress 0.01277127 
Run 38 stress 0.01177522 
Run 39 stress 0.01460151 
Run 40 stress 0.0124688 
Run 41 stress 0.01074402 
Run 42 stress 0.0107642 
Run 43 stress 0.009464615 
... New best solution
... Procrustes: rmse 0.00710886  max resid 0.02577806 
Run 44 stress 0.01433142 
Run 45 stress 0.01617776 
Run 46 stress 0.009562942 
... Procrustes: rmse 0.007442498  max resid 0.02525807 
Run 47 stress 0.01440775 
Run 48 stress 0.01442018 
Run 49 stress 0.01242664 
Run 50 stress 0.00946344 
... New best solution
... Procrustes: rmse 0.0002374678  max resid 0.0007088021 
... Similar to previous best
*** Best solution repeated 1 times

4.5.2 2. Do the antibiotics work?

4.5.2.1 Acclimation vs antibiotics

treat <- meta  %>% #antibiotics samples giving error when plotting
  filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

treat.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(treat))]
identical(sort(colnames(treat.counts)),sort(as.character(rownames(treat))))

treat_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")

4.5.2.2 Number of samples used

[1] 53
beta_div_richness_treat<-hillpair(data=treat.counts, q=0)
beta_div_neutral_treat<-hillpair(data=treat.counts, q=1)
beta_div_phylo_treat<-hillpair(data=treat.counts, q=1, tree=genome_tree)
beta_div_func_treat<-hillpair(data=treat.counts, q=1, dist=dist)
4.5.2.2.1 Richness
Df SumOfSqs R2 F Pr(>F)
time_point 1 2.033529 0.0957379 6.818959 0.001
species 1 2.327685 0.1095866 7.805342 0.001
individual 26 9.722175 0.4577167 1.253885 0.001
Residual 24 7.157208 0.3369589 NA NA
Total 52 21.240598 1.0000000 NA NA
4.5.2.2.2 Neutral
Df SumOfSqs R2 F Pr(>F)
time_point 1 2.168854 0.1067496 9.686182 0.001
species 1 3.090152 0.1520953 13.800733 0.001
individual 26 9.684304 0.4766554 1.663479 0.001
Residual 24 5.373891 0.2644996 NA NA
Total 52 20.317201 1.0000000 NA NA
4.5.2.2.3 Phylogenetic
Df SumOfSqs R2 F Pr(>F)
time_point 1 2.0671823 0.2208963 21.358337 0.001
species 1 0.8731487 0.0933035 9.021460 0.001
individual 26 4.0949687 0.4375828 1.627293 0.006
Residual 24 2.3228577 0.2482174 NA NA
Total 52 9.3581574 1.0000000 NA NA
4.5.2.2.4 Functional
Df SumOfSqs R2 F Pr(>F)
time_point 1 2.1253769 0.4009905 40.2325710 0.001
species 1 0.0070451 0.0013292 0.1333606 0.762
individual 26 1.9000410 0.3584768 1.3833480 0.185
Residual 24 1.2678545 0.2392035 NA NA
Total 52 5.3003175 1.0000000 NA NA
beta_richness_nmds_treat <- beta_div_richness_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = c("sample" = "Tube_code"))


beta_neutral_nmds_treat <- beta_div_neutral_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_treat <- beta_div_phylo_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_treat <- beta_div_func_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = join_by(sample == Tube_code))

4.5.3 3. Do the FMT work?

4.5.3.1 Comparison between FMT2 vs Post-FMT2

transplant3 <- meta  %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")

transplant3.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(transplant3))]
identical(sort(colnames(transplant3.counts)),sort(as.character(rownames(transplant3))))

transplant3_nmds <- sample_metadata %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")

4.5.3.2 Number of samples used

[1] 53
beta_div_richness_transplant3<-hillpair(data=transplant3.counts, q=0)
beta_div_neutral_transplant3<-hillpair(data=transplant3.counts, q=1)
beta_div_phylo_transplant3<-hillpair(data=transplant3.counts, q=1, tree=genome_tree)
beta_div_func_transplant3<-hillpair(data=transplant3.counts, q=1, dist=dist)
4.5.3.2.1 Richness
Df SumOfSqs R2 F Pr(>F)
species 1 0.9809525 0.0778442 5.775712 0.001
time_point 1 0.7092107 0.0562799 4.175734 0.001
type 1 1.4479860 0.1149060 8.525540 0.001
individual 25 6.9157236 0.5488022 1.628753 0.001
Residual 15 2.5476146 0.2021678 NA NA
Total 43 12.6014875 1.0000000 NA NA
4.5.3.2.2 Neutral
Df SumOfSqs R2 F Pr(>F)
species 1 1.1101711 0.0895457 7.879642 0.001
time_point 1 0.7363935 0.0593971 5.226687 0.001
type 1 1.9027421 0.1534740 13.505059 0.001
individual 25 6.5351386 0.5271203 1.855374 0.001
Residual 15 2.1133659 0.1704628 NA NA
Total 43 12.3978112 1.0000000 NA NA
4.5.3.2.3 Phylogenetic
Df SumOfSqs R2 F Pr(>F)
species 1 0.1341431 0.1005265 6.921553 0.001
time_point 1 0.1159136 0.0868654 5.980943 0.002
type 1 0.1423573 0.1066822 7.345390 0.001
individual 25 0.6512838 0.4880704 1.344205 0.070
Residual 15 0.2907075 0.2178554 NA NA
Total 43 1.3344053 1.0000000 NA NA
4.5.3.2.4 Functional
Df SumOfSqs R2 F Pr(>F)
species 1 -0.0031096 -0.0030953 -0.1177247 0.821
time_point 1 -0.0045864 -0.0045653 -0.1736359 0.811
type 1 0.0775554 0.0771987 2.9361453 0.138
individual 25 0.5385511 0.5360739 0.8155530 0.619
Residual 15 0.3962105 0.3943880 NA NA
Total 43 1.0046211 1.0000000 NA NA
beta_richness_nmds_transplant3 <- beta_div_richness_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_transplant3 <- beta_div_neutral_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == Tube_code))

beta_phylo_nmds_transplant3 <- beta_div_phylo_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_transplant3 <- beta_div_func_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == Tube_code))
p0<-beta_richness_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="top")

p1<-beta_neutral_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylo_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_func_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

4.5.4 4. Are there differences between the control and the treatment group?

4.5.4.1 After 1 week –> Post-FMT1

post1 <- meta %>%
  filter(time_point == "5_Post-FMT1")

post1.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post1))]
identical(sort(colnames(post1.counts)),sort(as.character(rownames(post1))))

post1_nmds <- sample_metadata %>%
  filter(time_point == "5_Post-FMT1")

4.5.4.2 Number of samples used

[1] 27
beta_div_richness_post1<-hillpair(data=post1.counts, q=0)
beta_div_neutral_post1<-hillpair(data=post1.counts, q=1)
beta_div_phylo_post1<-hillpair(data=post1.counts, q=1, tree=genome_tree)
beta_div_func_post1<-hillpair(data=post1.counts, q=1, dist=dist)
4.5.4.2.1 Richness
Df SumOfSqs R2 F Pr(>F)
species 1 0.6488906 0.0757032 2.115638 0.002
type 1 0.5615418 0.0655126 1.830847 0.009
Residual 24 7.3610760 0.8587842 NA NA
Total 26 8.5715084 1.0000000 NA NA
4.5.4.2.2 Neutral
Df SumOfSqs R2 F Pr(>F)
species 1 0.7984743 0.104299 3.065171 0.001
type 1 0.6051778 0.079050 2.323148 0.006
Residual 24 6.2519772 0.816651 NA NA
Total 26 7.6556293 1.000000 NA NA
4.5.4.2.3 Phylogenetic
Df SumOfSqs R2 F Pr(>F)
species 1 0.0700374 0.0635900 1.6594454 0.149
type 1 0.0184254 0.0167292 0.4365646 0.803
Residual 24 1.0129278 0.9196808 NA NA
Total 26 1.1013906 1.0000000 NA NA
4.5.4.2.4 Functional
Df SumOfSqs R2 F Pr(>F)
species 1 -0.0048931 -0.0058792 -0.1628739 0.837
type 1 0.1161466 0.1395538 3.8660898 0.082
Residual 24 0.7210172 0.8663254 NA NA
Total 26 0.8322707 1.0000000 NA NA
p0<-beta_richness_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="top")

p1<-beta_neutral_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")
grid.arrange(arrangeGrob(p0,p1,p2,p3))

4.5.4.3 After 2 weeks –>Post-FMT2

post2 <- meta %>%
  filter(time_point == "6_Post-FMT2")

post2.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post2))]
identical(sort(colnames(post2.counts)),sort(as.character(rownames(post2))))

post2_nmds <- sample_metadata %>%
  filter(time_point == "6_Post-FMT2")

4.5.4.4 Number of samples used

[1] 28
beta_div_richness_post2<-hillpair(data=post2.counts, q=0)
beta_div_neutral_post2<-hillpair(data=post2.counts, q=1)
beta_div_phylo_post2<-hillpair(data=post2.counts, q=1, tree=genome_tree)
beta_div_func_post2<-hillpair(data=post2.counts, q=1, dist=dist)
4.5.4.4.1 Richness
Df SumOfSqs R2 F Pr(>F)
species 1 0.8966107 0.1132233 3.515588 0.001
type 1 0.6463814 0.0816245 2.534445 0.002
Residual 25 6.3759661 0.8051521 NA NA
Total 27 7.9189582 1.0000000 NA NA
4.5.4.4.2 Neutral
Df SumOfSqs R2 F Pr(>F)
species 1 0.9398968 0.1224370 4.112304 0.001
type 1 1.0227481 0.1332297 4.474800 0.001
Residual 25 5.7139316 0.7443333 NA NA
Total 27 7.6765766 1.0000000 NA NA
4.5.4.4.3 Phylogenetic
Df SumOfSqs R2 F Pr(>F)
species 1 0.1200820 0.1454923 4.647187 0.001
type 1 0.0592745 0.0718175 2.293931 0.040
Residual 25 0.6459931 0.7826902 NA NA
Total 27 0.8253496 1.0000000 NA NA
4.5.4.4.4 Functional
Df SumOfSqs R2 F Pr(>F)
species 1 0.0116738 0.0168395 0.4229082 0.545
type 1 -0.0085273 -0.0123007 -0.3089208 0.893
Residual 25 0.6900904 0.9954612 NA NA
Total 27 0.6932368 1.0000000 NA NA
p0<-beta_richness_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="top")

p1<-beta_neutral_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")
grid.arrange(arrangeGrob(p0,p1,p2,p3))